Quantification of blood flow index in diffuse correlation spectroscopy using a robust deep learning method

Abstract. Significance Diffuse correlation spectroscopy (DCS) is a powerful, noninvasive optical technique for measuring blood flow. Traditionally the blood flow index (BFi) is derived through nonlinear least-square fitting the measured intensity autocorrelation function (ACF). However, the fitting process is computationally intensive, susceptible to measurement noise, and easily influenced by optical properties (absorption coefficient μa and reduced scattering coefficient μs′) and scalp and skull thicknesses. Aim We aim to develop a data-driven method that enables rapid and robust analysis of multiple-scattered light’s temporal ACFs. Moreover, the proposed method can be applied to a range of source–detector distances instead of being limited to a specific source–detector distance. Approach We present a deep learning architecture with one-dimensional convolution neural networks, called DCS neural network (DCS-NET), for BFi and coherent factor (β) estimation. This DCS-NET was performed using simulated DCS data based on a three-layer brain model. We quantified the impact from physiologically relevant optical property variations, layer thicknesses, realistic noise levels, and multiple source–detector distances (5, 10, 15, 20, 25, and 30 mm) on BFi and β estimations among DCS-NET, semi-infinite, and three-layer fitting models. Results DCS-NET shows a much faster analysis speed, around 17,000-fold and 32-fold faster than the traditional three-layer and semi-infinite models, respectively. It offers higher intrinsic sensitivity to deep tissues compared with fitting methods. DCS-NET shows excellent anti-noise features and is less sensitive to variations of μa and μs′ at a source–detector separation of 30 mm. Also, we have demonstrated that relative BFi (rBFi) can be extracted by DCS-NET with a much lower error of 8.35%. By contrast, the semi-infinite and three-layer fitting models result in significant errors in rBFi of 43.76% and 19.66%, respectively. Conclusions DCS-NET can robustly quantify blood flow measurements at considerable source–detector distances, corresponding to much deeper biological tissues. It has excellent potential for hardware implementation, promising continuous real-time blood flow measurements.


Introduction
Cerebral blood flow (CBF) is essential for monitoring metabolic oxygenation, 1,2 neurovascular coupling, 3,4 and metabolic response to functional stimuli. 5,6For example, CBF abnormalities are caused by ischemic strokes, 7 head trauma, 8 or brain injury. 9,10There are several blood flow measurement techniques, including computed tomography, 11 magnetic resonance imaging, 12 and positron emission tomography. 13However, although they are well-established, they cannot provide continuous, long-term measurements at the bedside.Laser Doppler flowmetry can measure microvascular blood flow but only probe shallow tissues. 14Doppler ultrasound techniques can only measure blood flow in larger vasculatures and are unsuitable for longitudinal monitoring for unstable probe orientations. 15Near-infrared diffuse optical methods are becoming popular in blood flow measurements as they are noninvasive, nonionized, portable, and faster.Among them is diffuse correlation spectroscopy (DCS), 16,17 using a laser with a long coherence length (>5 m) to illuminate tissue surfaces and collect remitted scattered light at a distance, typically 1 to 3 cm, away from the incident position.The scattered light from flowing red blood cells causes a speckle pattern fluctuating at a rate proportional to the flow rate.This blood-flow-dependent information can be quantified based on the normalized temporal intensity autocorrelation function (ACF) g 2 ðτÞ ≡ hIðtÞIðtþτÞi hIðtÞi 2 , where IðtÞ is the measured scattered light and τ is the correlation lag time. 17,18DCS can measure blood flow in vivo in small animals, 19,20 human brains, 16 and muscles. 21Traditionally, to derive blood flow index (BFi), the measured g 2 ðτÞ is fitted with a homogenous semi-infinite one-layer analytical model 22 or the Monte Carlo model. 235][26] However, treating biological tissues with a homogenous semi-infinite model is not quite realistic, as significant signal contamination from superficial tissue layers (e.g., scalp/skull) occurs when measuring deep flow in the brain.Research has been conducted to minimize the discrepancy, with the diffusion equation for layered geometries developed for fitting methods, including two- 27,28 and three-layer analytical models. 29,30Unfortunately, multilayer models highly rely on a priori knowledge of each layer's optical properties (namely the absorption coefficient μ a and reduced scattering coefficient μ 0 s ) and thickness to estimate blood flow within each layer.Commonly, layer optical properties and thicknesses are assumed from literature, and the errors in these assumed values can lead to significant errors in brain blood flow estimations.Additionally, the multilayer model is susceptible to measurement noise, especially for the three-layer model, although its accuracy in BFi estimations has been validated. 24,31Moreover, these methods are iterative and time-consuming.To overcome these limitations, the N'th-order linear (NL) algorithm, 32,33 least-absolute minimization (L1 norm), and the support vector regression (SVR) 34 were proposed.However, under the NL framework, BFi extraction is significantly influenced by the linear regression approach adopted. 34Although L1 norm and SVR are new approaches to processing DCS data, they are sensitive to signal deviations. 35,36Additionally, the BFi computing time is 28.07 and 52.93 s (using L1 norm and SVR, respectively), still slow for practical applications, particularly for real-time monitoring. 34eep learning, an increasingly popular method, has been widely applied to biomedical time sequence data, including electroencephalogram (EEG) and electrocardiogram (ECG), 37,38 but has yet to be broadly used in DCS.Very recently, Zhang et al. 39 proposed the first recurrent neural network (RNN) regression model to DCS, followed by 2D convolution neural networks (2D CNNs), 40 long short-term memory (LSTM), 41 and ConvGRU. 42LSTM, as a typical RNN structure, has proven stable and robust for quantifying relative blood flow in previous studies in phantom and in vivo experiments. 412D CNN, on the other hand, tends to require large training datasets for complex structures, demanding massive memory resources.ConvGRU, the newest deep learning method introduced to DCS, has also exhibited excellent performances in BFi extraction.
Nevertheless, all existing algorithms are designed for a single source-detector distance (ρ), corresponding to a specific depth in biological tissues.To accommodate a wider range of ρ, retraining the model becomes necessary.Inspired by a recently published one-dimensional convolutional neural network (1D CNN) 43 for fluorescence lifetime imaging (FLIM), we proposed the DCS neural network (DCS-NET) based on 1D CNN for quantifying the coherent factor β and BFi.
The primary objective of this work is to present and evaluate an artificial intelligence (AI) framework, called DCS-NET, in β and BFi estimations.We established the Monte Carlo simulation model based on the open-source tool Monte Carlo eXtreme (MCX) developed by Fang and Boas 44 to generate g 2 ðτÞ emulating experiment data.The DCS-NET training, validation, and testing datasets are from the semi-infinite geometry model. 22We investigated DCS-NET's performance on absolute BFi and relative BFI (rBFi)'s estimations and compared them with semiinfinite and three-layer model fitting methods.To best link our work with actual outcomes expected in practice, we modeled DCS measurement noise based on realistic experimental conditions, considering various noise levels controlled by the integration time (T int ).We define a metric that accounts for the intrinsic sensitivity of the brain blood flow and evaluate it between DCS-NET and traditional fitting methods.We also show BFi estimation errors induced by the inaccurate assumptions about layer optical properties and thicknesses when using fitting methods based on the semi-infinite and three-layer solutions of the correlation diffusion equation.Figure 1 summarizes the main concept of our work.All essential parameters are defined in Table 6 in the Appendix to facilitate our discussion.

DCS Theory
The transport of the unnormalized electric field auto-correlation function, G 1 ðρ; τÞ ≡ hE Ã ðρ; tÞ • Eðρ; tÞi, is well described by the correlation diffusion equation: 17,45 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 7 ; where k 0 ¼ 2πn∕λ is the wavenumber of light, n and λ are the refractive index and wavelength in the scattering medium, respectively.α is the fraction of dynamic photon scattering events in the medium.hΔr 2 ðτÞi is the mean squared displacement of scatterers in the turbid medium during a time interval τ.SðρÞ is the point source located at ρ; ρ is the source-detector distance.μ a and μ 0 s are the tissue's absorption and reduced scattering coefficients, respectively.For a semi-infinite medium, the solution of Eq. ( 1) using the extrapolated boundary condition for continuous-wave DCS is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 4 ; 5 5 0 , and z b ¼ 5∕ð3μ 0 s Þ to be consistent with Ref. 46.Previous studies have shown that the scatters' Brownian diffusion motion model 18,47 aligns well with in vivo DCS experiments, and therefore, the mean-squared displacement can be derived as hΔr 2 ðτÞi ¼ 6D b τ, where D b represents the effective diffusion coefficient.BFi in DCS is typically defined as αD b . 21,48 2 ðτÞ is linked to the normalized electric field auto-correlation function as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 4 ; 4 4 3 where β is a constant accounting for the collection setup, such as the number of detected speckles and the numerical aperture of the detection fiber.However, realistic biological tissues 49 show multiple layers with different physiological and optical properties.Using DCS to conduct in vivo CBF measurements, light must propagate through different layers, including the scalp and skull. 50,51Thus, layered analytical models have been proposed for BFi extraction.These include the two- 27,28 and three-layer analytical models. 24,30,31,52This study considers the three-layer analytical model, where a turbid medium consisting of N slabs is considered, as shown in Fig. 2(c).Each slab has its thickness, Δ p ¼ L p − L p−1 , p ¼ 1, 2, 3, where L 0;1;2;3 are the coordinates along the z-axis and μ a1;2;3 , and μ 0 s1;2;3 are absorption and scattering coefficients.To solve Eq. (1) in the layered medium (along z direction), we can use the Fourier transform Gðr; τÞ for the transverse coordinate ρ as where q is the radial spatial frequency.Equation (1) can then be rewritten as where Θ 2 ðpÞ ðq; τÞ ¼ 3μ aðpÞ μ 0 sðpÞ þ 6k 2 0 μ 02 sðpÞ D bðpÞ τ þ q 2 , z 0 ¼ 1∕μ 0 s1 , and p ¼ 1; 2; 3. We divided the top layer into two sublayers: layer 0 (0 < z < z 0 ) identified by p ¼ 0, and layer 1 (z 0 < z < Δ 1 ).Then, the solution of Eq. ( 5) inside the p'th layer (p ¼ 1; 2; 3) can be written as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 4 ; 1 3 5 Ĝðq; z; τÞ ¼ A ðpÞ expðΘ ðpÞ zÞ þ B ðpÞ expð−Θ ðpÞ zÞ; (6)   where A ðpÞ and B ðpÞ are coefficients for each layer determined by the boundary conditions where z 0 ∼ 1∕μ 0 s1 and z 3 ∼ 1∕μ 0 s3 are the extrapolation lengths accounting for internal reflections at the tissue surface (z ¼ 0) and the back surface (z ¼ L 3 ), respectively.D p ¼ c∕3μ 0 sðpÞ is the photon diffusion coefficient in layer p, and c is the speed of light.
Substituting Eq. ( 6) into Eq.( 7), A ðpÞ and B ðpÞ can be determined (p ¼ 1, 2, 3), and we obtain the solution of Eq. ( 5) at z ¼ 0 as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 7 ; 5 2 9 Ĝðq; z ¼ 0; τÞ ¼ Num Denom ; where Num and Denom (when p ¼ 3 and ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 7 ; 4 8 4 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 7 ; 4 3 0 Therefore, by performing the inverse Fourier transform of Eq. ( 8) with respect to q, the field ACF at z ¼ 0 can be written as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 1 7 ; 3 6 3 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 1 1 7 ; 3 0 1 where J 0 denotes the zero-order Bessel function of the first kind.The integral bound for q in Eq. ( 11) should theoretically be from 0 to þ∞.However, in practice, the numerical integration is performed with a limited range as ½0 mm −1 ; 30 mm −1 advised in Ref. 29.

Noise Models
This study evaluates the impact from noise on BFi and β.We employed a broadly accepted noise model proposed by Zhou et al. 53 The standard deviation (σðτÞ, noise) of g 2 ðτÞ is given as where T b is the bin width of the correlator, m is the bin index corresponding to τ. hni ≡ IT b is the average number of photons detected within the bin time, where I is the detected photon count rate, and T int is the integration time (e.g., measurement duration).Γ is the decay rate of g 2 ðτÞ, which is obtained from fitting the measured g 2 ðτÞ to the theoretical g 2 ðτÞ ≈ 1 þ β expð−ΓτÞ.Gaussian noise 54,55 was added to g 2 ðτÞ based on a statistical noise model to determine the noise (σðτÞ).Considering realistic photon budgets, the photon count rate at 785 nm was assumed to be 8.05 kcps. 55Three different noise levels were defined according to T int (= 1, 10, or 30 s).

Intrinsic Sensitivity Estimation
To evaluate the sensitivity to changes in blood flow in the deeper layer, we fixed the effective diffusion coefficient w is an integer and w ¼ 1;2; : : : ; 11.The physiological and optical parameters listed in Table 1 are taken as baseline conditions.Similar to Ref. 54, the intrinsic sensitivity (η H ) is defined as where BFi H and BFi 0 represent the estimated BFi (H ¼ D, S, or T, meaning DCS-NET, the semiinfinite, and three-layer fitting methods) for the perturbed and baseline conditions, respectively, and CBF perturb and CBF 0 are D b in layer 3 for the perturbed and baseline conditions, respectively.

Monte Carlo Simulations
We utilized a simplified model comprising three layers to emulate the scalp (5 mm), skull (7 mm), and brain (50 mm, large enough so that we can treat the medium as semi-infinite), respectively. 57All layers were assumed homogeneous, as demonstrated in Fig. 2(a), and their corresponding optical properties are summarized in Table 1.
MCX utilized an anisotropic factor (g) of 0.89 and a refractive index (n) of 1.37 44 for all layers.We launched 2 × 10 9 photons from a source with a diameter of 1 mm and set the detector radii to 0.13, 0.28, 0.45, 0.7, 1, and 1.5 mm for ρ ¼ 5, 10, 15, 20, 25, and 30 mm, respectively, recording data from multiple distances simultaneously.An example of the source and the detector was arranged as shown in Fig. 2(a).MCX records the path lengths and momentum transfer from the detected photons for obtaining the electric field ACF G 1 ðτÞ: where N p is the number of detected photons, N t is the number of tissue types (3 for our simulations), and Y s;i and L s;i stand for the total momentum transfer and the total path length of photon s in layer i, respectively.μ a;i is the absorption coefficient, and hΔr 2 ðτÞi i is the mean square displacement of the scattered particles in layer i.Here, hΔr 2 ðτÞi i ¼ 6D i τ, where D i is the effective diffusion coefficient of layer i.The simulated G 1 ðτÞ is normalized to G 1 ð0Þ, and then we can obtain g 2 ðτÞ using the Siegert relationship with β ¼ 0.5.In this simulation, the delay time 1 μs ≤ τ < 10;000 μs (127 data points) was used for g 2 ðτÞ.

Deep Learning Architecture Design
The structure of DCS-NET is shown in Fig. 3(a).DCS-NET takes g 2 ðτÞ to estimate β and BFi independently.DCS-NET consists of (1) a shared branch for temporal feature extraction and (2) two subsequent independent branches for estimating β and BFi, with a similar structure to the shared branch.The two CNN layers in the shared branch have a wider sliding window with a larger kernel size of 13 and a giant stride of 5.They are expected to capture more general features of the auto-correlation decay curves.The batch normalization (BN) layer 58 is employed after each convolutional layer.It reduces the shift of internal covariance and accelerates network training when processing normalized data.To implement feature pooling and effectively reconstruct β and BFi, we use a pointwise convolution layer with a kernel size of 1 after the convolutional neural network, followed by the activation function, the Sigmoid function.The model input is measured (here, we used data from MCX) g 2 ðτÞ, of which the size is 1 × 127.Both the estimated β and BFi have a size of 1 × 1.Note that the simulated g 2 ðτÞ was normalized to (0, 1] before being fed into the model.

Training Dataset Preparation
The training datasets can be easily obtained using synthetic data based on the homogenous semiinfinite analytical model, as shown in Fig. 2 We used an early stopping callback with 20 patient epochs to prevent overfitting.To match the realistic experiments, in the dataset, we set μ a ∈ Uð0.01; 1 mm −1 , μ 0 s ∈ Uð0.5; 1.6 mm −1 , β ∈ Uð0; 1, BFi ∈ ½10 −8 ; 10 −5 mm 2 ∕s, and ρ ∈ U½5; 30 mm, where U stands for a uniform distribution.g 2 ðτÞ training datasets contain noisy and noiseless (the noise model has been described in Ref. 53) ACFs, as shown in Fig. 3(c).The green, yellow, and red lines represent noisy g 2 ðτÞ, and the blue line represents noiseless g 2 ðτÞ.We used the optimizer Adam 59 for the training process, with the learning rate fixed at 1 × 10 −5 in the standard back-propagation.We used the mean square error loss function for updating the network by controlling the following problem: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 1 1 7 ; 1 3 4 where X is the network output (estimated BFi or β), and Y is the corresponding label (true BFi or β) in the i'th training pairs.F is the mapping function, ℘ is the trainable weights of our networks, and M is the number of training pairs.Figure 3 3 Results

Absolute BFi Recovery Versus Detection Depths
To investigate how the absolute BFi and β behave in terms of ρ among DCS-NET, semi-infinite, and three-layer fitting approaches, we generated g 2 ðτÞ via MCX Monte Carlo simulations for ρ ¼ 5, 10, 15, 20, 25, and 30 mm, as described in Sec.2.4.Table 1 shows all the relevant parameters used in MCX simulations.The absolute BFi in this study corresponds to the Brownian diffusion coefficient D b (assumed α ¼ 1).When using DCS-NET, g 2 ðτÞ was fed into the pre-trained model.For the semi-infinite fitting procedure, g 2 ðτÞ was fitted to Eqs. ( 2) and ( 3), and we assumed μ a ¼ 0.019 mm −1 , μ 0 s ¼ 1.099 mm −1 , for the brain layer (layer 3), as provided in Table 1.
We also fitted the simulated g 2 ðτÞ with the three-layer model, Eqs. ( 11) and ( 12), and where N τ is the number of sampled g 2 ðτÞ, and g 2 ðτÞ H is from Eq. (3) or Eq. ( 12).Fitting was performed on τ from 1 to 10;000 μs.Table 2 presents the true β and BFi and estimated β and BFi using DCS-NET, semi-infinite, and three-layer fitting methods.All input parameters for fitting are assumed as described above, and β GT ¼ 0.5.We define BFi D , BFi S , and BFi T (also β D , β S , and β T ) for DCS-NET, the semiinfinite, and three-layer fitting methods, respectively.We define ε BFi;D ð%Þ ¼ jBFi D − BFi GT j∕ BFi GT × 100%, where ε BFi;D is the BFi error with DCS-NET.Similarly, ε BFi;S and ε BFi;T are the BFi estimated errors with the semi-infinite and three-layer fitting methods.
Table 2 shows when the semi-infinite model is used, the estimated BFi is closer to layer 1 (αD b ¼ 1 × 10 −6 mm 2 ∕s), even for ρ ¼ 30 mm, suggesting that a homogenous fitting procedure is more sensitive to the superficial layers' dynamic properties.This finding is consistent with the results reported by Gagnon et al. 27 Using the three-layer fitting model, we obtained BFi T ¼ 7.15 × 10 −7 mm 2 ∕s, close to 1 × 10 −6 mm 2 ∕s when ρ ¼ 5 mm.This is because the mean light penetration depth is ∼ρ∕3 to ρ∕2. 19When ρ is small, most detected photons predominantly travel through layer 1.As ρ increases (ρ ≥ 10 mm), the estimated BFi decreases, reaching 5.63 × 10 −6 mm 2 ∕s at ρ ¼ 25 mm, with ε BFi;T of 6.17%.This is because as ρ increases, the detected photons penetrate inside the skull layer (αD b ¼ 0 mm 2 ∕s), resulting in an increased contribution of layer 2. This phenomenon is expected, because the three-layer modeling can remove the contribution from superficial layers 52 to obtain accurate BFi.Interestingly, when using DCS-NET, the estimated BFi increases as ρ increases, reaching 5.71 × 10 −6 mm 2 ∕s with ε BFi;D of 4.83% at ρ ¼ 30 mm.These results suggest that the AI model is capable of recognizing the depth.Regarding β estimation, there is no significant difference among the three methods.

Absolute BFi Recovery with Noise
Figure 3(c) displays the semi-infinite analytical example g 2 ðτÞ curves with noise using the model proposed by Zhou et al. 53 The curves were obtained with ρ ¼ 30 mm at different noise levels (T int ¼ 1; 10; 30 s), μ a ¼ 0.019 mm −1 , and μ 0 s ¼ 1.099 mm −1 with an assumed BFi ¼ 2 × 10 −7 mm 2 ∕s.To assess DCS-NET's performance in practical scenarios, we modified the Monte Carlo code to generate g 2 curves including noise according to Zhou et al.'s noise model. 53We generated 100 g 2 sets for each noise level (including noiseless).Still, we minimized Eq. ( 17) using the Levenberg-Marquardt optimization routine.We performed the residual analysis to assess the efficiency of the semi-infinite and three-layer models.We define the residual δ and resnorm (the squared 2-norm of the residual) ϵ as where q is the lag time index, and Q is the length of the time trace.fðβ; BFi; τ q Þ is the fitted g 2 ðτÞ obtained from fitting methods based on analytical models at the lag time τ q , and the corresponding true value is g 2 ðτ q Þ from MCX.The fitting results using the semi-infinite and threelayer analytical models are presented in Fig. 4, in which noisy g 2 ðτÞ curves from MCX (blue star-shaped) and fitted g 2 ðτÞ curves (red lines) at different noise levels are shown.We also calculated the mean BFi and β over 100 trials.As for β, we arrive at the same conclusion as Sec.3.1 that all three methods exhibit similar behaviors at the same noise level.A high noise level (T int ¼ 1 s) leads to a significant standard deviation, as shown in Fig. 5(a).Figure 5(b) shows the estimated BFi.The estimated BFi for the semi-infinite model deviates significantly from the ground truth.When using the three-layer fitting method, ε BFi;T is 82.30% at the lower noise level (T int ¼ 30 s).As the noise level increases, ε BFi;T also increases, with ε BFi;T reaching 390.10% at the high noise level (T int ¼ 1 s).Furthermore, a high noise level leads to a more significant standard deviation, indicating that BFi estimation is highly sensitive to noise when the three-layer fitting method is applied, in accordance with previous findings. 52In contrast, ε BFi;D (using DCS-NET) at a high noise level (T int ¼ 1 s) is 12.87%, whereas at a low noise level (T int ¼ 30 s), it is only 1.93%, indicating that DCS-NET is not susceptible to noise.Figure 5(b) also shows that when the three-layer fitting method is used, the BFi precision can be enhanced through increasing T int .
To compare the accuracy of the three different methods in quantifying rBFi, we defined the error in rBFi as ε rBFi;H ¼ jrBFi H − rBFi GT j∕rBFi GT × 100% (H ¼ D, S, or T), meaning the rBFi estimation error using DCS-NET, the semi-infinite, and three-layer fitting methods, respectively.We can observe that rBFi D (red star) is close to the true rBFi (blue solid line) with ε rBFi;D ranging from 0.15% to 8.35%.By contrast, the semi-infinite and three-layer methods result in more significant errors of 3.41% ≤ ε rBFi;S ≤ 43.76% and 0.36% ≤ ε rBFi;T ≤ 19.66%, respectively.As expected, the semi-infinite homogenous solution resulted in significant errors in rBFi, in agreement with Ref. 33.

Intrinsic Sensitivity
As described in Sec.2.3, the input D b in layer 3, denoted as CBF 0 ¼ 6 × 10 −6 mm 2 ∕s, serves as the base point, and its corresponding recovered BFi is denoted as BFi 0 .Similarly, we assigned αD b ¼ ½1 þ 0.05 × ðw − 1Þ × 6 × 10 −6 mm 2 ∕s (w is an integer; w ¼ 1;2; : : : ; 21), and it is referred to as the perturbed blood flow CBF perturb .We also define a perturbation level ζ ¼ ðCBF perturb − CBF 0 Þ∕CBF 0 × 100%.We calculated the corresponding BFi for αD b , and then used Eq. ( 14) to obtain η D , η S and η T .We considered physiological noise by utilizing the noise model described in Sec.2.2. Figure 7(a) shows the noiseless intrinsic sensitivity, demonstrating that DCS-NET exhibits η D > 71.34%.The intrinsic sensitivity reaches 2.5 × when ζ ¼ 20%, then decreases with ζ increasing.In comparison, the three-layer fitting method achieved η T ¼ 61.96%, whereas the semi-infinite fitting method yielded η S of only 14.12% on noiseless data.Figures 7(b)-7(d) illustrate sensitivity curves at various noise levels.Especially noteworthy are the instances where η D > 0 at T int ¼ 10 s and T int ¼ 30 s. Conversely, with the semi-infinite and three-layer fitting models, η predominantly assumes negative values, underscoring the considerable impact of measurement noise on sensitivity.Furthermore, the impact of measurement noise on the sensitivity overgrows, particularly for the three-layer fitting method, as apparent in Fig. 7(d).

BFi Extraction with Varied Optical Properties and Scalp/Skull Thicknesses
In practical applications, a patient's head parameters can vary significantly, and the ideal scenario is to measure them before conducting DCS measurements.However, it is not always straightforward, and we usually assume average values.However, we must evaluate the impact of assumed errors on BFi estimation.Since μ a and μ 0 s are typically unknown and have to be measured separately or taken from literature.We examined how μ a and μ 0 s of layer 3 (brain) impact BFi extraction.Changing the scalp/skull thickness also varies BFi, which can be observed using the multi-layered model fitting method.Here, we use the three-layer fitting method, and all BFi were obtained at ρ ¼ 30 mm.Additional details are presented in Table 3.

μ a variation
To study how μ a impacts BFi, we set μ a ¼ 0.011, 0.015, 0.019, 0.023, and 0.027 and μ 0 s ¼ 1.099 mm −1 in MCX.The baseline is at μ a ¼ 0.019 mm −1 , with AE20% and AE40% variation.In this case, two BFi groups were calculated.The first group was calculated assuming a constant μ a ¼ 0.019 mm −1 (0%), defined as μ a;m , and the calculated BFi is defined as BFi m .The second group was calculated using the known μ a set in MCX, which we considered as true μ a , and the corresponding calculated BFi is considered as BFi GT .

μ 0 s variation
Similarly, we conducted simulations with μ 0 s ¼ 0.666, 0.888, 1.110, 1.332, and 1.554 mm −1 and a fixed μ a ¼ 0.019 mm −1 to investigate how μ 0 s impacts BFi estimation.We define the estimated BFi as BFi m when μ 0 s ¼ 1.099 mm −1 (at 0%, defined as μ 0 s;m ).Additionally, BFi GT was calculated using the known μ 0 s set in MCX, considered as true μ 0 s .The mean and standard deviation of the estimated BFi (versus μ a ) over 100 trials are shown in Fig. 8(a).We also compare BFi m and BFi GT .The blue (BFi GT ) and green (BFi m ) dashed lines are for the semi-infinite model, whereas the red (BFi GT ) and purple (BFi m ) dashed line are for the three-layer model.The red solid (BFi GT ) and black dashed lines are for DCS-NET.Similarly, the BFi's mean and standard deviation (versus μ 0 s ) over 100 trials are shown in Fig. 8(b).
Figure 9 shows the BFi variation (in %) versus the μ a and μ 0 s variations (in %).The percentage error for μ a is defined as E μ a ¼ ½ μ a;m −μ a μ a × 100%.Similarly, we define the percentage error for μ 0 s as 8 and 9 show that E BFi is positively related to E μ a and negatively related to E μ 0 s for semi-infinite and three-layer fitting models, in good agreement with previous findings. 26,31On the other hand, E BFi curves obtained from DCS-NET are close and are not sensitive to E μ a and E μ 0 s .This result is expected, as from Eq. ( 2), μ 0 s should yield a more pronounced impact compared to μ a , primarily due to the second-order contribution from μ 0 s and μ 0 s ≫ μ a observed in biological tissues.Extreme E BFi examples are shown in Fig. 9, namely, a more extensive E μ a ∼ þ62% results in E BFi ∼ þ25% and E μ a ∼ −30% results in E BFi ∼ −10%.When E μ 0 s reaches +62%, E BFi reaches ∼ − 50% and E μ 0 s ∼ −30% gives E BFi ∼ þ70%.The results from the three-layer fitting model show similar behaviors.Namely, E BFi is positively related to E μ a and negatively related to E μ 0 s in layer 3, this result aligns well with the conclusions from Zhao et al.' conclusion. 31In contrast, DCS-NET only shows −1% ∼ þ5% in E BFi caused by E μ a and E μ 0 s (blue solid and brown solid lines for μ a and μ 0 s , respectively in Fig. 9), indicating that the variations in μ a and μ 0 s have negligible impact on BFi estimation.
Figure 10(a) presents BFi's mean value (represented by bar plots) and standard deviation (depicted by error bars) over 100 trials versus Δ 1 .The rightmost bar group represents the results obtained with Δ 1 ¼ 5 mm. Figure 10(b) shows BFi's mean value and standard deviation versus Δ 2 , the rightmost bar group represents the results obtained with Δ 2 ¼ 7 mm.Still, we can see that the semi-infinite model cannot provide accurate BFi at a deeper layer.When Δ 1 changed, ε BFi;D falls into 1.17% ∼ 8.33% [the bar group 1 in Fig. 10(a)] when using DCS-NET, whereas ε BFi;T falls into 4.30% ∼ 14.66% [the bar group 3 in Fig. 10(a)] using the three-layer fitting model, slightly larger than that using DCS-NET.However, ε BFi;T increases to 11.67% ∼ 16.05% when Δ 1 estimation error occurs using the three-layer fitting method [shown in the rightmost bar group in Fig. 10(a)].Whereas for the variation in Δ 2 , ε BFi;D falls into 0.33% ∼ 10.33% when DCS-NET is used [the bar group 1 in Fig. 10(b)], whereas ε BFi;T falls into 1.50% ∼ 13.33% when the three-layer fitting method is used [the bar group 3 in Fig. 10(b)].Both present similar accuracy.However, when Δ 2 is not accurate, ε BFi;T becomes more pronounced and reaches 41.09% ∼ 193.40% [the rightmost bar group in Fig. 10(b)].
Figure 11 shows the BFi variation (in %) versus the Δ 1 and Δ 2 variations (in %).The percentage error for Δ 1 is defined as Similarly, we define the percentage error for The BFi error (in %) caused by assumed error in E Δ 1 and E Δ 2 is defined as As it is commonly known, E Δ 1 and E Δ 2 cause a significant E BFi .Figures 10(a) and 10(b) demonstrate a positive correlation between E BFi and E Δ 1 (and E Δ 2 ).Furthermore, as observed in Fig. 11, E BFi resulting from E Δ 2 ranges from −176.41% to +43.68%.In contrast, E BFi caused by E Δ 1 ranges from −44.29% to +53.47%.This error range is significantly narrower than that caused by the skull thickness, agreeing with the findings in Ref. 31.For DCS-NET, E BFi caused by both Δ 1 and Δ 2 falls within the limited range of −6% to þ8%.

BFi Estimation Time
In addition, the BFi estimation time is also an important parameter, especially in real-time measurements, and Table 4 compares the three extraction methods.We record it for single decays and batch decays (e.g., 100 trials) at different noise levels.It is clear that DCS-NET is promising for real-time applications.All data reported in Table 4 are standard deviations and means for repeating three times after discarding the first few runs that usually take longer.The analysis were performed using the workstation (CPU: Intel(R) Core(TM) i9-10900X @3.70 GHz; Memory: 128 GB; graphics processing unit (GPU): NVIDIA Quadro RTX 5000).Our study shows that DCS-NET can robustly quantify DCS-based blood flow measurements.We used DCS-NET to analyze the ACFs generated from MCX.The proposed network is based on 1DCNN, 43 which is straightforward, quicker to train, and faster than high-dimension CNNs for time sequence analysis, such as FLIM data. 43,60To evaluate DCS-NET, we compared it with the semi-infinite, three-layer fitting methods by changing tissue optical properties (μ a and μ 0 s ), depths (related to ρ), and scalp/skull thicknesses (Δ 1 and Δ 2 ).BFi estimated by DCS-NET shows a small error range −1% ∼ þ5% induced by μ a and μ 0 s (see Fig. 9) and a slightly wider error range −6% ∼ þ8% induced by Δ 1 and Δ 2 (see Fig. 11).For rBFi, the error from DCS-NET (8.35%) is much less than that of the semi-infinite and three-layer fitting methods (43.76% and 19.66%, respectively).Moreover, DCS-NET yields more than 71.34% sensitivity to brain blood flow, whereas the semi-infinite and three-layer fitting methods yield 14.12% and 61.96%, respectively [Fig.7(a)].We considered measurement noise using a stochastic noise model 53 to reflect experimental realities.With DCS-NET, ε BFi;D is 12.87% at a high noise level (T int ¼ 1 s), whereas it increases to 390.10% when using the three-layer fitting method.At a low noise level (T int ¼ 30 s), the three-layer fitting model yields ε BFi;T of 82.30%, much worse than 1.93% obtained by DCS-NET, suggesting that DCS-NET is less sensitive to noise [see Fig. 5(b)].Figures 10(a) and 10(b) show that the three-layer analytical method (modeling the head, i.e., scalp, skull, and brain) can minimize the influence of extracerebral layers on measured DCS signals.However, this model requires a priori knowledge of layer optical properties and thicknesses.Therefore, accurately estimating scalp and skull thicknesses is required for reliable CBF estimation when using a three-layer analytical model.
Besides accuracy and robustness, the computational cost is a critical factor that impacts practical applications, especially for real-time monitoring.Table 4 reveals that it took 0.004 s for DCS-NET to quantify 100 g 2 curves with 127 data points.In contrast, it took 0.160 and 181.697 s, respectively, for the semi-infinite fitting and three-layer fitting procedures.For quantifying a single autocorrelation decay curve, it only took 0.001 s for DCS-NET.In contrast, it took 0.032 and 17.496 s, respectively, for the semi-infinite fitting and three-layer fitting procedures.DCS-NET is the fastest among the three, around 17,000-fold faster than the three-layer model and 32-fold faster than the semi-infinite model.
Table 5 lists existing deep learning methods applied to DCS techniques.It shows that DCS-NET's training is much faster than 2DCNNs, 40 approximately 140-fold faster.Although the remaining models, RNN, 39 LSTM, 41 and ConvGRU, 42 have fewer total layers, they are limited to a specific ρ.
Although DCS-NET is more robust than the semi-infinite and three-layer fitting methods, our study has several limitations.First, DCS-NET's training datasets were generated using the semi-infinite diffusion model as advised in Ref. 40.Nevertheless, this model does not consider scalp and skull thicknesses, which could potentially explain why the error range ð−6% ∼ þ8%Þ caused by Δ 1 and Δ 2 is much broader than that ð−1% ∼ þ5%Þ caused by μ a and μ 0 s (Figs. 9 and 11).The complexity of including training datasets generated from a layered model is beyond the scope of this study, given this report's already long length.In future, we will train new networks using datasets generated from a layered model, and alternatively, obtaining training datasets from in vivo measurements, as demonstrated in Refs.41 and 42 will also be considered.Second, current rBFi calculations do not consider variations in optical properties between the baseline and activation states.Indeed, μ a and μ 0 s in the brain can vary according to interventions (e.g., functional activation), which are recognized to impact perfusion.Failing to account for these changes could introduce additional uncertainties in rBFi measurements.Third, we did not include a comparison with the two-layered analytical model in this report; it may be worth further investigation.Fourth, as we all know, analytical fitting methods suffer from partial volume effects and recover only a fraction of the actual change; still, the relationship between the recovered change and the actual change remains linear.However, from Fig. 7(a), we can see the BFi values from DCS-NET reflect various degrees of the relative ground truth change according to the relative change; thus, they have a non-linear relationship with actual brain blood flow.This suggests processing data with our DCS-NET could result in non-physiological distortions.We will further investigate this and improve our network models in future studies.Finally, our study was solely conducted using simulation data.In the future, we will perform phantom and in vivo experiments to validate our findings.

Conclusion
We compared the proposed DCS-NET against the semi-infinite and the three-layer models for estimating β, BFi and rBFi.We used Monte Carlo simulations to validate their performances.This study evaluated the cerebral sensitivity using a deep learning method and the influence of scalp/skull thickness and μ a ∕μ 0 s variations on BFi extraction.Additionally, we examined the impact of noise.Our findings revealed that the homogenous model is sensitive to superficial layers.In contrast, the three-layer model performs better in estimating BFi in deeper layers but is more susceptible to measurement noise.
Furthermore, DCS-NET outperforms the semi-infinite and three-layer fitting models in rBFi recovery.Using DCS-NET, variations in μ a and μ 0 s have less impact on BFi, unlike variations in scalp and skull thicknesses, which show a more significant error in BFi.Moreover, iterative fitting methods are much slower and unsuitable for real-time "online" processing.In contrast, our DCS-NET is 32-fold faster than the semi-infinite model and 17,000-fold faster than the threelayer model, showing great potential for continuous real-time clinical applications.

Appendix
Table 6 shows all essential parameters used in the throughout article, ensuring accessibility to comprehensive details for interested readers.

Fig. 2
Fig. 2 Simulation layered model and analytical models.(a) A large slab representing a human brain consisting of three layers of the scalp (5 mm), skull (7 mm), and brain (50 mm).(b) The homogenous semi-infinite analytical model used for fitting methods and generating deep-learning training datasets.(c)Three-layer geometric scheme including the position of the source and detector, each layer has its own thickness Δ 1;2;3 and characterized by the absorption coefficient μ a1;2;3 and reduced scattering coefficient μ 0 s1;2;3 .
(b).Thus, according to Eqs. (2) and (3), 200,000 training datasets (200;000 × 127) were generated and split into the training (80%) and the validation (20%) groups.Each dataset consists of the input, g 2 ðτÞ, and its corresponding labels are BFi and β, which are the output.The training batch size is 128, with 800 training epochs.

Fig. 3
Fig. 3 Design and evaluation of the convolution neural network (CNN).(a) The proposed DCS-NET includes a CNN, BN, and sigmoid activation layers.The convolution layer parameters are the filter number × the kernel size × the stride.(b) Training and validation losses of DCS-NET.(c) g 2 ðτÞ with noise-free (blue), and with realistic noise added, assuming an 8.05 kcps at 785 nm at different noise levels with T int ¼ 1, 10, and 30 s.
(b) shows that the training and validation losses decrease rapidly and reach the plateau after 85 epochs.The training process's best score reaches a small value of 0.000725, indicating that the network is well trained as the estimated β and BFi are close to the ground truth.The model was conducted in Python using Pytorch with Intel (R) Core (TM) i9-10900KF CPU @3.70 GHz.

½g 2
and Δ 2 ¼ 7 mm.Meanwhile, we set β ¼ 0.3 and D b3 ¼ 2 × 10 −7 mm 2 ∕s as the initial guesses.For the fitting, we used NLSM (lsqcurvefitð•Þ in MATLAB with the Levenberg-Marquardt optimization) to minimize the unweighted least squares objective function, E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 ðτÞ MCX − g 2 ðτÞ H 2 ; H ¼ ðS; TÞ;

Figures 4 (
a)(i)-4(a)(iv) show the MCX-generated and fitted g 2 using the semi-infinite model, and they exhibit an increasing trend in δ, ranging from ð−0.0025; 0.0025Þ to ð−0.5; 0.5Þ, indicating that the semi-infinite method becomes inaccurate when the noise level increases.Additionally, ϵ reaches 3.02 when T int ¼ 1 s.Similar behaviors are observed in the three-layer fitting, as shown in Figs.4(b)(i)-4(b)(iv).

Fig. 4
Fig. 4 MCX-generated (scattered stars) and fitted (red solid lines) g 2 curves using semi-infinite and three-layer fitting methods.[(a) (i)-(iv), respectively] noisy MCX simulated data (scattered star-shaped) at different noise levels fitted with the semi-infinite homogeneous model; [(b) (i)-(iv)] noisy MCX-generated data fitted with the for the three-layer fitting procedure.The corresponding residual δ and resnorm ϵ curves are also included.

Fig. 5
Fig. 5 Estimated β by DCS-NET, semi-infinite, and three-layer fitting methods at different noise levels (T int ¼ 1, 10, and 30 s).The bar height means the average value for estimated BFi or β, the error bar means the standard deviations σ.(b) The estimated BFi by the three methods at different noise levels.The red dot line stands for the ground truth.(All the average values were obtained over 100 trials.)

Fig. 7
Fig. 7 (a) Intrinsic sensitivity on noiseless data.(b)-(d)The sensitivities for noise with T int ¼ 30 s, T int ¼ 10 s, T int ¼ 1 s, respectively.η is the intrinsic sensitivity that defined in Eq. (14), and ζ is the perturbation level in layer 3 (brain).Red, blue, and dark lines present η D , η T , and η S , respectively.The perturbation levels in the graphs start at ζ ¼ 20%.

Fig. 8
Fig.8(a) Estimated BFi versus μ a , the green and purple dashed lines are for BFi m assuming μ a ¼ 0.019 mm −1 , the red solid and black dashed lines are for BFi GT and BFi D , respectively, and the red and blue dashed lines are for BFi GT using the three-layer and semi-infinite fitting methods.(b) Estimated BFi versus μ 0 s , the green and purple dashed lines are for BFi m assuming μ 0 s ¼ 1.10 mm −1 , the red solid and black dashed lines are for BFi GT and BFi D , respectively, and the red and blue dashed lines are for BFi GT using the three-layer and semi-infinite fitting methods.

Fig. 10 (
Fig. 10 (a) BFi's mean value and standard deviation versus Δ 1 , and the rightmost bar group represents the results obtained with Δ 1 ¼ 5 mm.(b) BFi's mean value and standard deviation versus Δ 2 , and the rightmost bar group represents the results obtained with Δ 2 ¼ 7 mm.Each bar in the plot represents the average BFi over 100 trials calculated using three different methods, whereas the error bar stands for the standard deviation of BFi over 100 trials.

Table 1
56ysiological and optical parameters56at 785 nm in the human head model.

Table 2
BFi in the brain estimated using DCS-NET, homogeneous semi-infinite and three-layer fitting models.

Table 4
The BFi estimation time (with Matlab parfor for semi-infinite and three-layer fitting models).

Table 5
Comparison of existing AI methods for BFi estimation.
Notes: the training parameters of RNN and CNN(2D) are not given in the literature; we calculate them according to the structure shown in the literature.

Table 6
Essential parameters list.The percentage error of the assumed μ a Defined in Sec.3.5 s The percentage error of the assumed μ 0 s Defined in Sec.3.5 ρ Source-detector distance

Table 6 (
Continued).Estimated BFi when assumed μ a , μ 0 s , Δ 1 , Δ 2 are constant at 0% Defined in Sec.3.5 Blood flow index estimated by the semi-infinite fitting method Defined in Sec.3.1 BFi T Blood flow index estimated by the three-layer fitting method Defined in Sec.3.1 E BFi (%) The BFi error (in %) between BFi m and BFi GT Defined in Sec.3.5 ε BFi;D Error percentage of BFi using DCS-NET Defined in Sec.3.1 ε BFi;S Error percentage of BFi using the semi-infinite fitting method Defined in Sec.3.1 wInteger, w ¼ 1; 2; 3 δ Residual from fitting procedures Defined in Eq. (18)ϵ Resnorm (the square 2-norm of the residual) Defined in Eq.(18)